Interpretable scRNA-seq Trajectory DE with scLANE

Author
Affiliation

Jack Leary

University of Florida

Published

September 20, 2023

Introduction

In this tutorial we’ll walk through a basic trajectory differential expression analysis. We’ll use the scLANE R package, which we developed with the goal of providing accurate and biologically interpretable models of expression over pseudotime. At the end are a list of references we used in developing the method & writing the accompanying manuscript, as well as the poster I presented at ENAR 2023 in Nashville.

Libraries

If you haven’t already, install the development version (currently v0.7.3) of scLANE from the GitHub repository.

Code
remotes::install_github("jr-leary7/scLANE")

Next, we’ll load the packages we need to process, analyze, & visualize our data.

Code
library(dplyr)           # data manipulation
library(scLANE)          # trajectory DE 
library(Seurat)          # scRNA-seq tools
library(ggplot2)         # plot creation 
library(patchwork)       # plot combination
library(slingshot)       # pseudotime estimation
library(reticulate)      # Python interface
library(ComplexHeatmap)  # heatmaps
rename <- dplyr::rename

Helper Functions

We’ll also define a couple utilities to make our plots cleaner to read & easier to make.

Code
theme_umap <- function(base.size = 14) {
  ggplot2::theme(axis.ticks = ggplot2::element_blank(), 
                 axis.text = ggplot2::element_blank(), 
                 plot.subtitle = ggplot2::element_text(face = "italic", size = 11), 
                 plot.caption = ggplot2::element_text(face = "italic", size = 11))
}
guide_umap <- function(key.size = 4) {
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = key.size, alpha = 1)))
}

And consistent color palettes will make our plots easier to understand.

Code
palette_cluster <- paletteer::paletteer_d("ggsci::default_jama")
palette_celltype <- paletteer::paletteer_d("ggsci::category20_d3")
palette_heatmap <- paletteer::paletteer_d("wesanderson::Zissou1")

Data

We’ll load the well-known pancreatic endocrinogenesis data from Bastidas-Ponce et al (2019), which comes with the scVelo Python library & has been used in several pseudotime inference / RNA velocity method papers as a good benchmark dataset due to the simplicity of the underlying trajectory manifold.

Code
import scvelo as scv
adata = scv.datasets.pancreas()

  0%|          | 0.00/50.0M [00:00<?, ?B/s]
  1%|          | 400k/50.0M [00:00<00:13, 4.00MB/s]
  2%|1         | 0.98M/50.0M [00:00<00:09, 5.18MB/s]
  3%|3         | 1.56M/50.0M [00:00<00:09, 5.56MB/s]
  4%|4         | 2.12M/50.0M [00:00<00:08, 5.62MB/s]
  5%|5         | 2.66M/50.0M [00:00<00:08, 5.61MB/s]
  6%|6         | 3.12M/50.0M [00:00<00:09, 5.32MB/s]
  7%|7         | 3.59M/50.0M [00:00<00:09, 5.18MB/s]
  8%|8         | 4.05M/50.0M [00:00<00:09, 5.00MB/s]
  9%|9         | 4.54M/50.0M [00:00<00:09, 5.01MB/s]
 10%|9         | 4.96M/50.0M [00:01<00:09, 4.80MB/s]
 11%|#         | 5.43M/50.0M [00:01<00:09, 4.84MB/s]
 12%|#1        | 5.89M/50.0M [00:01<00:09, 4.77MB/s]
 13%|#2        | 6.34M/50.0M [00:01<00:09, 4.73MB/s]
 14%|#3        | 6.87M/50.0M [00:01<00:09, 4.95MB/s]
 15%|#4        | 7.43M/50.0M [00:01<00:08, 5.12MB/s]
 16%|#5        | 7.98M/50.0M [00:01<00:08, 5.28MB/s]
 17%|#7        | 8.52M/50.0M [00:01<00:08, 5.34MB/s]
 18%|#8        | 9.05M/50.0M [00:01<00:08, 5.36MB/s]
 19%|#9        | 9.62M/50.0M [00:01<00:07, 5.55MB/s]
 20%|##        | 10.2M/50.0M [00:02<00:07, 5.61MB/s]
 21%|##1       | 10.7M/50.0M [00:02<00:07, 5.46MB/s]
 22%|##2       | 11.2M/50.0M [00:02<00:07, 5.43MB/s]
 23%|##3       | 11.7M/50.0M [00:02<00:07, 5.30MB/s]
 24%|##4       | 12.2M/50.0M [00:02<00:07, 5.21MB/s]
 25%|##5       | 12.6M/50.0M [00:02<00:07, 5.03MB/s]
 26%|##6       | 13.1M/50.0M [00:02<00:07, 4.96MB/s]
 27%|##7       | 13.6M/50.0M [00:02<00:07, 5.03MB/s]
 28%|##8       | 14.2M/50.0M [00:02<00:07, 5.31MB/s]
 30%|##9       | 14.8M/50.0M [00:02<00:06, 5.40MB/s]
 31%|###       | 15.3M/50.0M [00:03<00:06, 5.50MB/s]
 32%|###1      | 15.9M/50.0M [00:03<00:06, 5.54MB/s]
 33%|###2      | 16.5M/50.0M [00:03<00:06, 5.62MB/s]
 34%|###4      | 17.0M/50.0M [00:03<00:06, 5.73MB/s]
 35%|###5      | 17.5M/50.0M [00:03<00:06, 5.58MB/s]
 36%|###5      | 18.0M/50.0M [00:03<00:06, 5.41MB/s]
 37%|###7      | 18.6M/50.0M [00:03<00:05, 5.51MB/s]
 38%|###8      | 19.1M/50.0M [00:03<00:05, 5.57MB/s]
 39%|###9      | 19.7M/50.0M [00:03<00:05, 5.58MB/s]
 40%|####      | 20.2M/50.0M [00:03<00:05, 5.66MB/s]
 41%|####1     | 20.8M/50.0M [00:04<00:05, 5.58MB/s]
 43%|####2     | 21.3M/50.0M [00:04<00:05, 5.71MB/s]
 44%|####3     | 21.8M/50.0M [00:04<00:06, 4.67MB/s]
 45%|####4     | 22.3M/50.0M [00:04<00:05, 4.88MB/s]
 46%|####5     | 22.9M/50.0M [00:04<00:05, 5.13MB/s]
 47%|####6     | 23.4M/50.0M [00:04<00:05, 5.30MB/s]
 48%|####7     | 23.9M/50.0M [00:04<00:05, 5.29MB/s]
 49%|####8     | 24.4M/50.0M [00:04<00:05, 5.27MB/s]
 50%|####9     | 24.9M/50.0M [00:04<00:05, 5.09MB/s]
 51%|#####     | 25.4M/50.0M [00:05<00:04, 5.18MB/s]
 52%|#####1    | 26.0M/50.0M [00:05<00:04, 5.32MB/s]
 53%|#####2    | 26.5M/50.0M [00:05<00:04, 5.42MB/s]
 54%|#####3    | 27.0M/50.0M [00:05<00:04, 5.27MB/s]
 55%|#####4    | 27.5M/50.0M [00:05<00:04, 5.16MB/s]
 56%|#####5    | 27.9M/50.0M [00:05<00:04, 5.07MB/s]
 57%|#####6    | 28.4M/50.0M [00:05<00:04, 4.95MB/s]
 58%|#####7    | 28.8M/50.0M [00:05<00:04, 4.83MB/s]
 59%|#####8    | 29.3M/50.0M [00:05<00:04, 4.81MB/s]
 59%|#####9    | 29.8M/50.0M [00:05<00:04, 4.79MB/s]
 60%|######    | 30.2M/50.0M [00:06<00:04, 4.85MB/s]
 61%|######1   | 30.7M/50.0M [00:06<00:04, 4.84MB/s]
 62%|######2   | 31.2M/50.0M [00:06<00:04, 4.84MB/s]
 63%|######3   | 31.6M/50.0M [00:06<00:04, 4.82MB/s]
 64%|######4   | 32.1M/50.0M [00:06<00:03, 4.79MB/s]
 65%|######5   | 32.6M/50.0M [00:06<00:03, 4.83MB/s]
 66%|######6   | 33.1M/50.0M [00:06<00:03, 4.83MB/s]
 67%|######6   | 33.5M/50.0M [00:06<00:03, 4.84MB/s]
 68%|######7   | 34.0M/50.0M [00:06<00:03, 4.83MB/s]
 69%|######8   | 34.5M/50.0M [00:06<00:03, 4.91MB/s]
 70%|######9   | 35.0M/50.0M [00:07<00:03, 5.05MB/s]
 71%|#######   | 35.5M/50.0M [00:07<00:02, 5.11MB/s]
 72%|#######1  | 36.0M/50.0M [00:07<00:02, 5.01MB/s]
 73%|#######2  | 36.4M/50.0M [00:07<00:02, 4.93MB/s]
 74%|#######3  | 36.9M/50.0M [00:07<00:02, 4.92MB/s]
 75%|#######4  | 37.4M/50.0M [00:07<00:02, 4.95MB/s]
 76%|#######5  | 37.9M/50.0M [00:07<00:02, 5.00MB/s]
 77%|#######6  | 38.5M/50.0M [00:07<00:02, 5.24MB/s]
 78%|#######7  | 39.0M/50.0M [00:07<00:02, 5.23MB/s]
 79%|#######8  | 39.4M/50.0M [00:08<00:02, 5.15MB/s]
 80%|#######9  | 40.0M/50.0M [00:08<00:01, 5.29MB/s]
 81%|########  | 40.5M/50.0M [00:08<00:01, 5.24MB/s]
 82%|########1 | 40.9M/50.0M [00:08<00:01, 5.02MB/s]
 83%|########2 | 41.5M/50.0M [00:08<00:01, 5.24MB/s]
 84%|########4 | 42.1M/50.0M [00:08<00:01, 5.45MB/s]
 85%|########5 | 42.6M/50.0M [00:08<00:01, 5.53MB/s]
 86%|########6 | 43.2M/50.0M [00:08<00:01, 5.56MB/s]
 87%|########7 | 43.7M/50.0M [00:08<00:01, 5.50MB/s]
 88%|########8 | 44.2M/50.0M [00:08<00:01, 5.57MB/s]
 90%|########9 | 44.8M/50.0M [00:09<00:00, 5.64MB/s]
 91%|######### | 45.4M/50.0M [00:09<00:00, 5.76MB/s]
 92%|#########1| 45.9M/50.0M [00:09<00:00, 5.71MB/s]
 93%|#########2| 46.5M/50.0M [00:09<00:00, 5.63MB/s]
 94%|#########3| 47.0M/50.0M [00:09<00:00, 5.52MB/s]
 95%|#########5| 47.6M/50.0M [00:09<00:00, 5.73MB/s]
 96%|#########6| 48.1M/50.0M [00:09<00:00, 5.72MB/s]
 97%|#########7| 48.7M/50.0M [00:09<00:00, 5.77MB/s]
 98%|#########8| 49.3M/50.0M [00:09<00:00, 5.90MB/s]
100%|#########9| 49.8M/50.0M [00:09<00:00, 5.79MB/s]
100%|##########| 50.0M/50.0M [00:09<00:00, 5.26MB/s]

The AnnData object contains data that we’ll need to extract, specifically the counts matrices (stored in AnnData.layers) and the cell-level metadata (which is in AnnData.obs).

Code
adata
AnnData object with n_obs × n_vars = 3696 × 27998
    obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
    var: 'highly_variable_genes'
    uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'distances', 'connectivities'

Conversion from Python

The reticulate package allows us to pass the counts matrices & metadata from Python back to R. We’ll use the spliced mRNA counts as our default assay, and also define a new assay containing the total (spliced + unspliced) mRNA in each cell. Lastly, we remove genes with non-zero spliced mRNA in 3 or fewer cells. Note: while downloading this dataset requires a Python installation as well as the installation of the scVelo Python library (and its dependencies), running scLANE is done purely in R & requires no Python whatsoever.

Code
spliced_counts <- Matrix::Matrix(t(as.matrix(py$adata$layers["spliced"])), sparse = TRUE)
unspliced_counts <- Matrix::Matrix(t(as.matrix(py$adata$layers["unspliced"])), sparse = TRUE)
rna_counts <- spliced_counts + unspliced_counts
colnames(rna_counts) <- colnames(spliced_counts) <- colnames(unspliced_counts) <- py$adata$obs_names$to_list()
rownames(rna_counts) <- rownames(spliced_counts) <- rownames(unspliced_counts) <- py$adata$var_names$to_list()
spliced_assay <- CreateAssayObject(counts = spliced_counts)
spliced_assay@key <- "spliced_"
unspliced_assay <- CreateAssayObject(counts = unspliced_counts)
unspliced_assay@key <- "unspliced_"
rna_assay <- CreateAssayObject(counts = rna_counts)
rna_assay@key <- "rna_"
meta_data <- py$adata$obs %>% 
             mutate(cell_name = rownames(.), .before = 1) %>% 
             rename(celltype = clusters, 
                    celltype_coarse = clusters_coarse) %>% 
             mutate(nCount_spliced = colSums(spliced_counts), 
                    nFeature_spliced = colSums(spliced_counts > 0), 
                    nCount_unspliced = colSums(unspliced_counts), 
                    nFeature_unspliced = colSums(unspliced_counts > 0), 
                    nCount_rna = colSums(rna_counts), 
                    nFeature_rna = colSums(rna_counts > 0))
seu <- CreateSeuratObject(counts = spliced_assay, 
                          assay = "spliced", 
                          project = "Mm_Panc_Endo", 
                          meta.data = meta_data)
seu@assays$unspliced <- unspliced_assay
seu@assays$RNA <- rna_assay
seu <- seu[rowSums(seu@assays$spliced) > 3, ]

Preprocessing

We preprocess the counts using a typical pipeline with QC, normalization & scaling, dimension reduction, and graph-based clustering via the Leiden algorithm.

Code
seu <- PercentageFeatureSet(seu, 
                            pattern = "^mt-", 
                            col.name = "percent_mito", 
                            assay = "spliced") %>% 
       PercentageFeatureSet(pattern = "^Rp[sl]", 
                            col.name = "percent_ribo", 
                            assay = "spliced") %>% 
       NormalizeData(assay = "spliced", verbose = FALSE) %>% 
       NormalizeData(assay = "unspliced", verbose = FALSE) %>% 
       NormalizeData(assay = "RNA", verbose = FALSE) %>% 
       FindVariableFeatures(assay = "spliced", 
                            nfeatures = 3000, 
                            verbose = FALSE) %>% 
       ScaleData(assay = "spliced", 
                 vars.to.regress = c("percent_mito", "percent_ribo"), 
                 model.use = "poisson", 
                 verbose = FALSE) %>% 
       RunPCA(assay = "spliced", 
              npcs = 30, 
              approx = TRUE, 
              seed.use = 312, 
              verbose = FALSE) %>% 
       RunUMAP(reduction = "pca", 
               dims = 1:30, 
               n.components = 2, 
               metric = "cosine", 
               seed.use = 312, 
               verbose = FALSE) %>% 
       FindNeighbors(reduction = "pca", 
                     k.param = 30,
                     nn.method = "annoy", 
                     annoy.metric = "cosine", 
                     verbose = FALSE) %>% 
       FindClusters(algorithm = 4, 
                    method = "igraph", 
                    resolution = 0.5, 
                    random.seed = 312, 
                    verbose = FALSE)

Let’s visualize the results on our UMAP embedding. The clustering generally agrees with the celltype labels, though there is some overclustering in the ductal cells & underclustering in the mature endocrine celltypes.

Code
p0 <- DimPlot(seu, 
              group.by = "seurat_clusters", 
              pt.size = 1, 
              cols = alpha(palette_cluster, 0.75)) + 
      labs(color = "Leiden Cluster") + 
      theme_scLANE() + 
      theme_umap() + 
      theme(plot.title = element_blank(), 
            axis.title = element_blank(), 
            axis.line.x = element_blank()) + 
      guide_umap()
p1 <- DimPlot(seu, 
              group.by = "celltype", 
              pt.size = 1, 
              cols = alpha(palette_celltype, 0.75)) + 
      labs(x = "UMAP 1", 
           y = "UMAP 2", 
           color = "Celltype") + 
      theme_scLANE() + 
      theme_umap() + 
      theme(plot.title = element_blank()) + 
      guide_umap()
p2 <- (p0 / p1) +
      plot_layout(guides = "collect") + 
      plot_annotation(title = "Murine pancreatic endocrinogenesis", 
                      theme = theme_scLANE())
p2

Trajectory Inference

Pseudotime Estimation

We’ll start by fitting a trajectory using the slingshot R package. We define cluster 5 as the starting cluster, since in this case we’re already aware of the dataset’s underlying biology. After generating the estimates for each cell, we rescale the ordering to be defined on \([0, 1]\). This has no effect on the trajectory DE results however, and is mostly an aesthetic choice.

Code
sling_res <- slingshot(Embeddings(seu, "umap"),
                       start.clus = "5",
                       clusterLabels = seu$seurat_clusters, 
                       approx_points = 500)
sling_pt <- slingPseudotime(sling_res) %>% 
            as.data.frame() %>% 
            magrittr::set_colnames(c("PT")) %>% 
            mutate(PT = (PT - min(PT)) / (max(PT) - min(PT)))
seu <- AddMetaData(seu, 
                   metadata = sling_pt, 
                   col.name = "sling_pt")

Let’s visualize the results on our UMAP embedding. They match what we would expect (knowing the biological background of the data), with ductal cells at the start of the process and endocrine celltypes such as alpha, beta, & delta cells at the end of it.

Code
p3 <- Embeddings(seu, "umap") %>% 
      as.data.frame() %>% 
      magrittr::set_colnames(c("UMAP_1", "UMAP_2")) %>% 
      mutate(PT = sling_pt$PT) %>% 
      ggplot(aes(x = UMAP_1, y = UMAP_2, color = PT)) + 
      geom_point(size = 1, alpha = 0.75) + 
      labs(color = "Pseudotime") + 
      scale_color_gradientn(colors = palette_heatmap, 
                            labels = scales::label_number(accuracy = 0.01)) + 
      theme_scLANE() + 
      theme_umap() + 
      theme(axis.title = element_blank(), 
            axis.line.x = element_blank())
p4 <- (p3 / p1) + 
      plot_layout(guides = "collect") + 
      plot_annotation(title = "Estimated cell ordering from Slingshot", 
                      theme = theme_scLANE())
p4

Trajectory Differential Expression

Next, we prepare the primary inputs to scLANE: our Seurat object with the spliced counts set as the default assay, a dataframe containing our estimated pseudotime ordering, a vector of size factors to use as an offset in each model, and a set of genes whose dynamics we want to model. scLANE parallelizes over genes in order to speed up the computation at the expense of using a little more memory. The models are fit using NB GLMs with optimal spline knots identified empirically, and differential expression is quantified using a likelihood ratio test of the fitted model vs. a constant (intercept-only) model. In practice, genes designated as HVGs are usually the best candidates for modeling, so we choose the top 3,000 HVGs as our input. Note: the testing of the HVG set on its own is also justified by the reality that almost all trajectories are inferred using some sort of dimension-reduced space, and those embeddings are nearly universally generated using a set of HVGs. As such, genes not included in the HVG set actually have no direct relationship with the estimated trajectory, & it’s generally safe to exclude them from trajectory analyses.

Code
top3k_hvg <- HVFInfo(seu) %>% 
             arrange(desc(variance.standardized)) %>% 
             slice_head(n = 3000) %>% 
             rownames(.)
cell_offset <- createCellOffset(seu)
scLANE_res <- testDynamic(seu, 
                          pt = sling_pt, 
                          genes = top3k_hvg, 
                          size.factor.offset = cell_offset, 
                          n.cores = 4, 
                          track.time = TRUE)
[1] "testDynamic evaluated 3000 genes across 1 lineage in 22.121 mins"

After tidying the TDE results with getResultsDE(), we pull a sample of 6 genes from the results & display their test statistics. By default, any gene with an adjusted p-value less than 0.01 is predicted to be dynamic, though this threshold can be easily adjusted.

Code
scLANE_res_tidy <- getResultsDE(scLANE_res)
select(scLANE_res_tidy, 
       Gene, 
       Test_Stat, 
       P_Val, 
       P_Val_Adj,
       Gene_Dynamic_Overall) %>% 
  mutate(Gene_Dynamic_Overall = if_else(Gene_Dynamic_Overall == 1, "Dynamic", "Static")) %>% 
  with_groups(Gene_Dynamic_Overall, 
              slice_sample, 
              n = 3) %>% 
  kableExtra::kbl(digits = 5, 
                  booktabs = TRUE, 
                  col.names = c("Gene", "LRT Statistic", "P-value", "Adj. P-value", "Predicted Gene Status")) %>% 
  kableExtra::kable_classic(full_width = FALSE, "hover")
Gene LRT Statistic P-value Adj. P-value Predicted Gene Status
Chmp2b 201.17610 0.00000 0.00000 Dynamic
B2m 288.15731 0.00000 0.00000 Dynamic
Timeless 281.09855 0.00000 0.00000 Dynamic
Abhd11os 18.63898 0.00009 0.03685 Static
Ctss 5.28626 0.07114 1.00000 Static
Klk10 -96.01652 1.00000 1.00000 Static

Next, we can use the plotModels() function to visualize the fitted models from scLANE and compare them to other modeling methods. The gene Neurog3 is strongly associated with epithelial cell differentiation, and indeed we see a very clear, nonlinear transcriptional dynamic across pseudotime for that gene. A traditional GLM fails to capture that nonlinearity, and while a GAM fits the trend smoothly, it fails to capture the sharpness of the transcriptional switch that occurs halfway through the trajectory. Only the scLANE model accurately models the rapid upregulation and equally swift downregulation of Neurog3 over pseudotime, in addition to identifying the cutpoint in pseudotime at which downregulation begins.

Code
p5 <- plotModels(scLANE_res, 
                 gene = "Neurog3", 
                 pt = sling_pt, 
                 expr.mat = seu, 
                 size.factor.offset = cell_offset, 
                 plot.null = FALSE) + 
      scale_color_manual(values = c("forestgreen"))
p5

Downstream analysis

Gene dynamics plots

Using the getFittedValues() function allows us to generate predictions from the models we fit, which we then use to visualize the dynamics of a few genes that are known to be strongly associated with the differentiation of immature cells into mature endocrine phenotypes. For all four genes, the fitted models show knots chosen in the area of pseudotime around the pre-endocrine cells. This tells us that these driver genes are being upregulated in precursor celltypes & are driving differentiation into the mature celltypes such as alpha & beta cells, after which the genes are downregulated.

Code
p6 <- getFittedValues(scLANE_res, 
                      genes = c("Chga", "Chgb", "Fev", "Cck"), 
                      pt = sling_pt, 
                      expr.mat = seu, 
                      size.factor.offset = cell_offset, 
                      cell.meta.data = select(seu@meta.data, celltype, celltype_coarse)) %>% 
      ggplot(aes(x = pt, y = rna_log1p)) + 
      facet_wrap(~gene, 
                 ncol = 2, 
                 scales = "free_y") + 
      geom_point(aes(color = celltype), size = 1, alpha = 0.75) + 
      geom_vline(data = data.frame(gene = "Chga", knot = unique(scLANE_res$Chga$Lineage_A$MARGE_Slope_Data$Breakpoint)), 
                 mapping = aes(xintercept = knot), 
                 linetype = "dashed", 
                 color = "grey20") + 
      geom_vline(data = data.frame(gene = "Chgb", knot = unique(scLANE_res$Chgb$Lineage_A$MARGE_Slope_Data$Breakpoint)), 
                 mapping = aes(xintercept = knot), 
                 linetype = "dashed", 
                 color = "grey20") + 
      geom_vline(data = data.frame(gene = "Cck", knot = unique(scLANE_res$Cck$Lineage_A$MARGE_Slope_Data$Breakpoint)), 
                 mapping = aes(xintercept = knot), 
                 linetype = "dashed", 
                 color = "grey20") + 
      geom_vline(data = data.frame(gene = "Fev", knot = unique(scLANE_res$Fev$Lineage_A$MARGE_Slope_Data$Breakpoint)), 
                 mapping = aes(xintercept = knot), 
                 linetype = "dashed", 
                 color = "grey20") + 
      geom_ribbon(aes(ymin = scLANE_ci_ll_log1p, ymax = scLANE_ci_ul_log1p), 
                  linewidth = 0, 
                  fill = "grey70", 
                  alpha = 0.9) + 
      geom_line(aes(y = scLANE_pred_log1p), 
                color = "black", 
                linewidth = 0.75) + 
      scale_x_continuous(labels = scales::label_number(accuracy = 0.01)) + 
      scale_color_manual(values = palette_celltype) + 
      labs(x = "Pseudotime", 
           y = "Normalized Expression", 
           title = "Endrocrinogenesis driver genes across pseudotime", 
           subtitle = "scLANE piecewise negative binomial GLMs") + 
      theme_scLANE() + 
      theme(legend.title = element_blank(), 
            strip.text.x = element_text(face = "italic"), 
            plot.subtitle = element_text(face = "italic", size = 11)) + 
      guide_umap()
p6

On the other hand, if we use additive models the “peak” of expression is placed among the mature endocrine celltypes - which doesn’t make biological sense if we know that these genes are driving that process of differentiation. This can of course be tweaked by changing the degree or degrees of freedom of the underlying basis spline, but choosing a “best” value for those hyperparameters can be difficult, whereas scLANE identifies optimal parameters internally by default. In addition, the knots chosen by scLANE for each gene can be informative with respect to the underlying biology, whereas the knots from GAMs are evenly spaced at quantiles & carry no biological interpretation.

Code
p7 <- getFittedValues(scLANE_res, 
                      genes = c("Chga", "Chgb", "Fev", "Cck"), 
                      pt = sling_pt, 
                      expr.mat = seu, 
                      size.factor.offset = cell_offset, 
                      cell.meta.data = select(seu@meta.data, celltype, celltype_coarse)) %>% 
      mutate(rna_raw = rna / size_factor, .before = 7) %>% 
      with_groups(gene, 
                  mutate, 
                  GAM_fitted_link = predict(nbGAM(expr = rna_raw, 
                                                  pt = sling_pt, 
                                                  Y.offset = cell_offset, 
                                                  spline.df = 3)), 
                  GAM_se_link = predict(nbGAM(expr = rna_raw, 
                                              pt = sling_pt, 
                                              Y.offset = cell_offset, 
                                              spline.df = 3), se.fit = T)[[2]]) %>% 
      mutate(GAM_pred = exp(GAM_fitted_link) * cell_offset, 
             GAM_ci_ll = exp(GAM_fitted_link - qnorm(0.975) * GAM_se_link) * cell_offset, 
             GAM_ci_ul = exp(GAM_fitted_link + qnorm(0.975) * GAM_se_link) * cell_offset, 
             GAM_pred_log1p = log1p(GAM_pred), 
             GAM_ci_ll_log1p = log1p(GAM_ci_ll), 
             GAM_ci_ul_log1p = log1p(GAM_ci_ul)) %>% 
      ggplot(aes(x = pt, y = rna_log1p)) + 
      facet_wrap(~gene, 
                 ncol = 2, 
                 scales = "free_y") + 
      geom_point(aes(color = celltype), size = 1, alpha = 0.75) + 
      geom_ribbon(aes(ymin = GAM_ci_ll_log1p, ymax = GAM_ci_ul_log1p), 
                  linewidth = 0, 
                  fill = "grey70", 
                  alpha = 0.9) + 
      geom_line(aes(y = GAM_pred_log1p), 
                color = "black", 
                linewidth = 0.75) + 
      scale_x_continuous(labels = scales::label_number(accuracy = 0.01)) + 
      scale_color_manual(values = palette_celltype) + 
      labs(x = "Pseudotime", 
           y = "Normalized Expression", 
           title = "Endrocrinogenesis driver genes across pseudotime", 
           subtitle = "Cubic basis spline negative binomial GAMs") + 
      theme_scLANE() + 
      theme(legend.title = element_blank(), 
            strip.text.x = element_text(face = "italic"), 
            plot.subtitle = element_text(face = "italic", size = 11)) + 
      guide_umap()
p7

Distribution of knot locations

Let’s take a broader view of the dataset by examining the distribution of adaptively chosen knots from our models. We limit the analysis to the set of genes determined to be dynamic.

Code
dyn_genes <- filter(scLANE_res_tidy, Gene_Dynamic_Overall == 1) %>% 
             pull(Gene)
knot_df <- purrr::imap(scLANE_res[dyn_genes], 
                       \(x, y) {
                         data.frame(
                           gene = y, 
                           knot = x$Lineage_A$MARGE_Slope_Data$Breakpoint
                         )
                       }) %>% 
           purrr::reduce(rbind)

We’ll plot a histogram of the knot values along with a ridgeplot of the pseudotime distribution for each celltype. We see that the majority of the selected knots are placed at the beginning of the trajectory, around where the ductal cells transition into endocrine progenitors. A smaller set of knots is placed about halfway through the trajectory, which we’ve annotated as the point at which pre-endocrine cells begin differentiating into mature endocrine phenotypes.

Code
p8 <- ggplot(knot_df, aes(x = knot)) + 
      geom_density(fill = "deepskyblue3", 
                   alpha = 0.75, 
                   color = "deepskyblue4", 
                   linewidth = 1) + 
      scale_x_continuous(limits = c(0, 1), labels = scales::label_number(accuracy = 0.1)) + 
      labs(x = "Knot Location") + 
      theme_scLANE() + 
      theme(axis.title.y = element_blank(), 
            axis.text.y = element_blank(), 
            axis.ticks.y = element_blank())
p9 <- data.frame(celltype = seu$celltype, 
                 pt = seu$sling_pt) %>% 
      ggplot(aes(x = pt, y = celltype, fill = celltype, color = celltype)) + 
      ggridges::geom_density_ridges(alpha = 0.75, size = 1, scale = 0.95) + 
      scale_x_continuous(labels = scales::label_number(accuracy = 0.1), limits = c(0, 1)) + 
      scale_fill_manual(values = palette_celltype) + 
      scale_color_manual(values = palette_celltype) + 
      labs(x = "Pseudotime") + 
      theme_scLANE() + 
      theme(axis.title.y = element_blank(), 
            legend.title = element_blank()) + 
      guide_umap()
p10 <- (p8 / p9) + 
       plot_layout(heights = c(1/4, 3/4)) + 
       plot_annotation(title = "Distribution of adaptively-chosen knots from scLANE", 
                       theme = theme_scLANE())
p10
Picking joint bandwidth of 0.0184

Dynamic gene clustering

We can extract a matrix of fitted values using smoothedCountsMatrix(); here we focus on the top 2,000 most dynamic genes, with the goal of identifying clusters of similarly-expressed genes. After reducing dimensionality with PCA, we cluster the genes using the Leiden algorithm & embed the genes in two dimensions with UMAP.

Code
smoothed_counts <- smoothedCountsMatrix(scLANE_res, 
                                        pt = sling_pt, 
                                        genes = dyn_genes[1:2000], 
                                        size.factor.offset = cell_offset, 
                                        n.cores = 2)
set.seed(312)
smoothed_counts_pca <- irlba::prcomp_irlba(t(smoothed_counts$Lineage_A), 
                                           n = 30, 
                                           center = TRUE, 
                                           scale. = TRUE)
smoothed_counts_umap <- uwot::umap(smoothed_counts_pca$x, 
                                   n_components = 2, 
                                   metric = "cosine", 
                                   n_neighbors = 20, 
                                   init = "spectral")
smoothed_counts_snn <- bluster::makeSNNGraph(smoothed_counts_pca$x, 
                                             k = 20, 
                                             type = "jaccard", 
                                             BNPARAM = BiocNeighbors::AnnoyParam(distance = "Cosine"))
smoothed_counts_clust <- igraph::cluster_leiden(smoothed_counts_snn, 
                                                objective_function = "modularity", 
                                                resolution_parameter = 0.3)
gene_clust_df <- data.frame(gene = colnames(smoothed_counts$Lineage_A), 
                            pc1 = smoothed_counts_pca$x[, 1], 
                            pc2 = smoothed_counts_pca$x[, 2], 
                            umap1 = smoothed_counts_umap[, 1], 
                            umap2 = smoothed_counts_umap[, 2], 
                            leiden = as.factor(smoothed_counts_clust$membership - 1L))

First we’ll visualize the gene clusters on the first two PCs.

Code
p11 <- ggplot(gene_clust_df, aes(x = pc1, y = pc2, color = leiden)) + 
       geom_point(size = 1, alpha = 0.75) + 
       labs(x = "PC 1", 
            y = "PC 2", 
            color = "Leiden Cluster", 
            title = "Unsupervised Clustering of Dynamic Genes", 
            subtitle = "Top 2,000 TDE genes after PCA") +
       paletteer::scale_color_paletteer_d("ggsci::default_igv") + 
       theme_scLANE() + 
       theme_umap() + 
       guide_umap()
p11

The UMAP embedding shows that even with the relatively small number of genes, clear patterns are visible.

Code
p12 <- ggplot(gene_clust_df, aes(x = umap1, y = umap2, color = leiden)) + 
       geom_point(size = 1, alpha = 0.75) + 
       labs(x = "UMAP 1", 
            y = "UMAP 2", 
            color = "Leiden Cluster", 
            title = "Unsupervised Clustering of Dynamic Genes", 
            subtitle = "Top 2,000 TDE genes after PCA") +
       paletteer::scale_color_paletteer_d("ggsci::default_igv") + 
       theme_scLANE() + 
       theme_umap() + 
       guide_umap()
p12

Expression cascades

We can also plot a heatmap of the dynamic genes; this requires a bit of setup, for which we’ll use the ComplexHeatmap package. We scale each gene, and clip values to be on \([-6, 6]\). The columns (cells) of the heatmap are ordered by estimated pseudotime, and the rows (genes) are ordered by expression peak.

Code
col_anno_df <- select(seu@meta.data, 
                      cell_name, 
                      celltype, 
                      sling_pt) %>% 
               mutate(celltype = as.factor(celltype)) %>% 
               arrange(sling_pt)
gene_order <- sortGenesHeatmap(smoothed_counts$Lineage_A, pt.vec = sling_pt$PT)
heatmap_mat <- t(scale(smoothed_counts$Lineage_A))
heatmap_mat[heatmap_mat > 6] <- 6
heatmap_mat[heatmap_mat < -6] <- -6
colnames(heatmap_mat) <- seu$cell_name
heatmap_mat <- heatmap_mat[, col_anno_df$cell_name]
heatmap_mat <- heatmap_mat[gene_order, ]
palette_celltype_hm <- as.character(palette_celltype[1:length(unique(seu$celltype))])
names(palette_celltype_hm) <- levels(col_anno_df$celltype)
col_anno <- HeatmapAnnotation(Celltype = col_anno_df$celltype, 
                              Pseudotime = col_anno_df$sling_pt, 
                              col = list(Celltype = palette_celltype_hm, 
                                         Pseudotime = circlize::colorRamp2(seq(0, 1, by = 0.25), palette_heatmap)),
                              show_legend = TRUE, 
                              show_annotation_name = FALSE, 
                              gap = unit(1, "mm"), 
                              border = TRUE)
palette_cluster_hm <- as.character(paletteer::paletteer_d("ggsci::default_igv")[1:length(unique(gene_clust_df$leiden))])
names(palette_cluster_hm) <- as.character(unique(gene_clust_df$leiden))
row_anno <- HeatmapAnnotation(Cluster = as.factor(gene_clust_df$leiden), 
                              col = list(Cluster = palette_cluster_hm), 
                              show_legend = TRUE, 
                              show_annotation_name = FALSE, 
                              annotation_legend_param = list(title = "Gene\nCluster"), 
                              gap = unit(1, "mm"), 
                              border = TRUE, 
                              which = "row")

The heatmap shows clear dynamic patterns across pseudotime; these patterns are often referred to as expression cascades, and represent periodic up- and down-regulation of different gene programs during the course of the underlying biological process.

Code
Heatmap(matrix = heatmap_mat, 
        name = "Spliced\nmRNA", 
        col = circlize::colorRamp2(colors = viridis::inferno(50), 
                                   breaks = seq(min(heatmap_mat), max(heatmap_mat), length.out = 50)), 
        cluster_columns = FALSE,
        width = 12, 
        height = 6, 
        column_title = "Dynamic genes across pseudotime in murine pancreatic endocrinogenesis",
        column_title_gp = gpar(fontface = "bold"), 
        cluster_rows = FALSE,
        top_annotation = col_anno, 
        left_annotation = row_anno, 
        show_column_names = FALSE, 
        show_row_names = FALSE, 
        use_raster = TRUE,
        raster_by_magick = TRUE, 
        raster_quality = 5)
Loading required namespace: magick

Enrichment analysis & gene programs

Using our gene clusters & the {gprofiler2} package, we run an enrichment analysis against the biological process (BP) set of gene ontologies.

Code
gene_clust_list <- purrr::map(unique(gene_clust_df$leiden), \(x) { 
  filter(gene_clust_df, leiden == x) %>% 
  pull(gene)
}) 
names(gene_clust_list) <- paste0("Leiden_", unique(gene_clust_df$leiden))
enrich_res <- gprofiler2::gost(gene_clust_list, 
                               organism = "mmusculus", 
                               ordered_query = FALSE, 
                               multi_query = FALSE, 
                               sources = "GO:BP", 
                               significant = TRUE)

A look at the top 3 most-significant GO terms for each gene cluster reveals heterogeneous functionalities across groups of genes. Cluster 1 in particular looks interesting; it shows significant enrichment of peptide regulation-related processes.

Code
mutate(enrich_res$result, 
       query = gsub("Leiden_", "", query)) %>% 
  rename(cluster = query) %>% 
  with_groups(cluster, 
              slice_head,
              n = 3) %>% 
  select(cluster, term_name, p_value, term_size, query_size, intersection_size, term_id) %>% 
  kableExtra::kbl(digits = 5, 
                  booktabs = TRUE, 
                  caption = "<i>Top 3 Biological Process GO Terms per Cluster<\\i>", 
                  col.names = c("Gene Cluster", "Term Name", "Adj. P-value", "Term Size", 
                                "Query Size", "Intersection Size", "Term ID")) %>% 
  kableExtra::kable_classic(c("hover"), full_width = FALSE)
Top 3 Biological Process GO Terms per Cluster
Gene Cluster Term Name Adj. P-value Term Size Query Size Intersection Size Term ID
0 peptide hormone secretion 0 313 253 28 GO:0030072
0 regulation of peptide hormone secretion 0 262 253 26 GO:0090276
0 peptide secretion 0 320 253 28 GO:0002790
1 organonitrogen compound metabolic process 0 6389 230 117 GO:1901564
1 negative regulation of biological process 0 6031 230 112 GO:0048519
1 regulation of apoptotic process 0 1656 230 53 GO:0042981
2 cell cycle 0 1802 212 143 GO:0007049
2 cell cycle process 0 1240 212 126 GO:0022402
2 mitotic cell cycle 0 870 212 103 GO:0000278
3 animal organ development 0 3278 277 88 GO:0048513
3 system development 0 4076 277 99 GO:0048731
3 anatomical structure development 0 6222 277 128 GO:0048856
4 secretion 0 1084 161 32 GO:0046903
4 secretion by cell 0 908 161 29 GO:0032940
4 export from cell 0 977 161 29 GO:0140352
5 behavior 0 765 199 30 GO:0007610
5 cell-cell signaling 0 1741 199 43 GO:0007267
5 synaptic signaling 0 928 199 30 GO:0099536
6 system development 0 4076 179 70 GO:0048731
6 multicellular organism development 0 4828 179 77 GO:0007275
6 anatomical structure development 0 6222 179 88 GO:0048856
7 multicellular organism development 0 4828 285 117 GO:0007275
7 anatomical structure development 0 6222 285 135 GO:0048856
7 developmental process 0 6808 285 141 GO:0032502

We isolate the genes from cluster 1, then create a per-cell module score for genes in that set. In effect, this will allow us to associate certain gene programs with certain celltypes.

Code
peptide_gene_program <- filter(gene_clust_df, leiden == 1) %>% 
                        pull(gene)
seu <- AddModuleScore(seu, 
                      features = list(peptide = peptide_gene_program), 
                      assay = "spliced", 
                      name = "peptide_program_score", 
                      seed = 312)

Visualizing the scores on our UMAP embedding shows us that the peptide program is most highly-enriched in mature endocrine cells. This makes sense biologically as mature endocrine celltypes’ primary roles are to produce peptides such as glucagon (alpha cells), insulin (beta cells), somatostatin (ductal cells), and pancreatic polypeptide (gamma cells).

Code
p13 <- Embeddings(seu, "umap") %>% 
       as.data.frame() %>% 
       magrittr::set_colnames(c("UMAP_1", "UMAP_2")) %>% 
       mutate(peptide_program_score = seu$peptide_program_score1) %>% 
       ggplot(aes(x = UMAP_1, y = UMAP_2, color = peptide_program_score)) + 
       geom_point(size = 1, alpha = 0.75) + 
       labs(color = "Gene Program Score") + 
       scale_color_gradientn(colors = palette_heatmap, 
                             labels = scales::label_number(accuracy = 0.01)) + 
       theme_scLANE() + 
       theme_umap() + 
       theme(axis.title = element_blank(), 
             axis.line.x = element_blank())
p14 <- (p13 / p1) + 
       plot_layout(guides = "collect") + 
       plot_annotation(title = "Enrichment of peptide regulation gene program", 
                       theme = theme_scLANE())
p14

Conclusions

Hopefully this vignette has been a useful introduction to running the scLANE software and using its outputs to help better understand biology at single-cell resolution. If you have questions about how the models work or are interpreted, software issues, or simply want to compare results feel free to open an issue on the GitHub repository, or reach out via email to .

Session Info

Code
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.1 (2022-06-23)
 os       macOS Big Sur ... 10.16
 system   x86_64, darwin17.0
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/New_York
 date     2023-09-20
 pandoc   2.19.2 @ /usr/local/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date (UTC) lib source
 abind                  1.4-5      2016-07-21 [1] CRAN (R 4.2.0)
 assertthat             0.2.1      2019-03-21 [1] CRAN (R 4.2.0)
 backports              1.4.1      2021-12-13 [1] CRAN (R 4.2.0)
 bigassertr             0.1.5      2021-07-08 [1] CRAN (R 4.2.0)
 bigparallelr           0.3.2      2021-10-02 [1] CRAN (R 4.2.0)
 bigstatsr              1.5.6      2022-02-03 [1] CRAN (R 4.2.0)
 Biobase              * 2.56.0     2022-04-26 [1] Bioconductor
 BiocGenerics         * 0.42.0     2022-04-26 [1] Bioconductor
 BiocNeighbors          1.14.0     2022-04-26 [1] Bioconductor
 BiocParallel           1.30.3     2022-06-05 [1] Bioconductor
 bitops                 1.0-7      2021-04-24 [1] CRAN (R 4.2.0)
 bluster                1.6.0      2022-04-26 [1] Bioconductor
 boot                   1.3-28     2021-05-03 [1] CRAN (R 4.2.1)
 broom                  1.0.0      2022-07-01 [1] CRAN (R 4.2.0)
 broom.mixed            0.2.9.4    2022-04-17 [1] CRAN (R 4.2.0)
 Cairo                  1.6-0      2022-07-05 [1] CRAN (R 4.2.0)
 circlize               0.4.15     2022-05-10 [1] CRAN (R 4.2.0)
 cli                    3.6.1      2023-03-23 [1] CRAN (R 4.2.0)
 clue                   0.3-61     2022-05-30 [1] CRAN (R 4.2.0)
 cluster                2.1.4      2022-08-22 [1] CRAN (R 4.2.0)
 coda                   0.19-4     2020-09-30 [1] CRAN (R 4.2.0)
 codetools              0.2-18     2020-11-04 [1] CRAN (R 4.2.1)
 colorspace             2.0-3      2022-02-21 [1] CRAN (R 4.2.0)
 ComplexHeatmap       * 2.12.1     2022-08-09 [1] Bioconductor
 cowplot                1.1.1      2020-12-30 [1] CRAN (R 4.2.0)
 crayon                 1.5.1      2022-03-26 [1] CRAN (R 4.2.0)
 data.table             1.14.2     2021-09-27 [1] CRAN (R 4.2.0)
 DBI                    1.1.3      2022-06-18 [1] CRAN (R 4.2.0)
 DelayedArray           0.22.0     2022-04-26 [1] Bioconductor
 DelayedMatrixStats     1.18.0     2022-04-26 [1] Bioconductor
 deldir                 1.0-6      2021-10-23 [1] CRAN (R 4.2.0)
 digest                 0.6.29     2021-12-01 [1] CRAN (R 4.2.0)
 doParallel             1.0.17     2022-02-07 [1] CRAN (R 4.2.0)
 dplyr                * 1.0.9      2022-04-28 [1] CRAN (R 4.2.0)
 ellipsis               0.3.2      2021-04-29 [1] CRAN (R 4.2.0)
 emmeans                1.8.0      2022-08-05 [1] CRAN (R 4.2.0)
 estimability           1.4.1      2022-08-05 [1] CRAN (R 4.2.0)
 evaluate               0.16       2022-08-09 [1] CRAN (R 4.2.0)
 fansi                  1.0.3      2022-03-24 [1] CRAN (R 4.2.0)
 farver                 2.1.1      2022-07-06 [1] CRAN (R 4.2.0)
 fastmap                1.1.0      2021-01-25 [1] CRAN (R 4.2.0)
 fitdistrplus           1.1-8      2022-03-10 [1] CRAN (R 4.2.0)
 flock                  0.7        2016-11-12 [1] CRAN (R 4.2.0)
 forcats                0.5.2      2022-08-19 [1] CRAN (R 4.2.0)
 foreach                1.5.2      2022-02-02 [1] CRAN (R 4.2.0)
 furrr                  0.3.1      2022-08-15 [1] CRAN (R 4.2.0)
 future                 1.27.0     2022-07-22 [1] CRAN (R 4.2.0)
 future.apply           1.9.0      2022-04-25 [1] CRAN (R 4.2.0)
 gamlss                 5.4-3      2022-04-24 [1] CRAN (R 4.2.0)
 gamlss.data            6.0-2      2021-11-07 [1] CRAN (R 4.2.0)
 gamlss.dist            6.0-5      2022-08-28 [1] CRAN (R 4.2.1)
 geeM                   0.10.1     2018-06-18 [1] CRAN (R 4.2.0)
 generics               0.1.3      2022-07-05 [1] CRAN (R 4.2.0)
 GenomeInfoDb         * 1.32.3     2022-08-09 [1] Bioconductor
 GenomeInfoDbData       1.2.8      2022-08-29 [1] Bioconductor
 GenomicRanges        * 1.48.0     2022-04-26 [1] Bioconductor
 GetoptLong             1.0.5      2020-12-15 [1] CRAN (R 4.2.0)
 ggh4x                  0.2.2      2022-08-14 [1] CRAN (R 4.2.0)
 ggplot2              * 3.4.2      2023-04-03 [1] CRAN (R 4.2.0)
 ggrepel                0.9.1      2021-01-15 [1] CRAN (R 4.2.0)
 ggridges               0.5.3      2021-01-08 [1] CRAN (R 4.2.0)
 glm2                 * 1.2.1      2018-08-11 [1] CRAN (R 4.2.0)
 glmmTMB                1.1.5      2022-11-16 [1] CRAN (R 4.2.0)
 GlobalOptions          0.1.2      2020-06-10 [1] CRAN (R 4.2.0)
 globals                0.16.1     2022-08-28 [1] CRAN (R 4.2.1)
 glue                   1.6.2      2022-02-24 [1] CRAN (R 4.2.0)
 goftest                1.2-3      2021-10-07 [1] CRAN (R 4.2.0)
 gprofiler2             0.2.1      2021-08-23 [1] CRAN (R 4.2.0)
 gridExtra              2.3        2017-09-09 [1] CRAN (R 4.2.0)
 gtable                 0.3.0      2019-03-25 [1] CRAN (R 4.2.0)
 here                   1.0.1      2020-12-13 [1] CRAN (R 4.2.0)
 highr                  0.9        2021-04-16 [1] CRAN (R 4.2.0)
 htmltools              0.5.3      2022-07-18 [1] CRAN (R 4.2.0)
 htmlwidgets            1.5.4      2021-09-08 [1] CRAN (R 4.2.0)
 httpuv                 1.6.5      2022-01-05 [1] CRAN (R 4.2.0)
 httr                   1.4.4      2022-08-17 [1] CRAN (R 4.2.0)
 ica                    1.0-3      2022-07-08 [1] CRAN (R 4.2.0)
 igraph                 1.3.4      2022-07-19 [1] CRAN (R 4.2.0)
 IRanges              * 2.30.1     2022-08-18 [1] Bioconductor
 irlba                  2.3.5      2021-12-06 [1] CRAN (R 4.2.0)
 iterators              1.0.14     2022-02-05 [1] CRAN (R 4.2.0)
 jsonlite               1.8.0      2022-02-22 [1] CRAN (R 4.2.0)
 kableExtra             1.3.4      2021-02-20 [1] CRAN (R 4.2.0)
 KernSmooth             2.23-20    2021-05-03 [1] CRAN (R 4.2.1)
 knitr                  1.40       2022-08-24 [1] CRAN (R 4.2.0)
 labeling               0.4.2      2020-10-20 [1] CRAN (R 4.2.0)
 later                  1.3.0      2021-08-18 [1] CRAN (R 4.2.0)
 lattice                0.20-45    2021-09-22 [1] CRAN (R 4.2.1)
 lazyeval               0.2.2      2019-03-15 [1] CRAN (R 4.2.0)
 leiden                 0.4.2      2022-05-09 [1] CRAN (R 4.2.0)
 lifecycle              1.0.3      2022-10-07 [1] CRAN (R 4.2.0)
 listenv                0.8.0      2019-12-05 [1] CRAN (R 4.2.0)
 lme4                   1.1-30     2022-07-08 [1] CRAN (R 4.2.0)
 lmtest                 0.9-40     2022-03-21 [1] CRAN (R 4.2.0)
 magick                 2.7.3      2021-08-18 [1] CRAN (R 4.2.0)
 magrittr             * 2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
 MASS                   7.3-58.1   2022-08-03 [1] CRAN (R 4.2.0)
 Matrix                 1.4-1      2022-03-23 [1] CRAN (R 4.2.1)
 MatrixGenerics       * 1.8.1      2022-06-30 [1] Bioconductor
 matrixStats          * 0.62.0     2022-04-19 [1] CRAN (R 4.2.0)
 mgcv                   1.8-40     2022-03-29 [1] CRAN (R 4.2.1)
 mime                   0.12       2021-09-28 [1] CRAN (R 4.2.0)
 miniUI                 0.1.1.1    2018-05-18 [1] CRAN (R 4.2.0)
 minqa                  1.2.4      2014-10-09 [1] CRAN (R 4.2.0)
 multcomp               1.4-20     2022-08-07 [1] CRAN (R 4.2.0)
 munsell                0.5.0      2018-06-12 [1] CRAN (R 4.2.0)
 mvtnorm                1.1-3      2021-10-08 [1] CRAN (R 4.2.0)
 nlme                   3.1-159    2022-08-09 [1] CRAN (R 4.2.0)
 nloptr                 2.0.3      2022-05-26 [1] CRAN (R 4.2.0)
 numDeriv               2016.8-1.1 2019-06-06 [1] CRAN (R 4.2.0)
 paletteer              1.5.0      2022-10-19 [1] CRAN (R 4.2.0)
 parallelly             1.32.1     2022-07-21 [1] CRAN (R 4.2.0)
 patchwork            * 1.1.2      2022-08-19 [1] CRAN (R 4.2.0)
 pbapply                1.5-0      2021-09-16 [1] CRAN (R 4.2.0)
 pillar                 1.8.1      2022-08-19 [1] CRAN (R 4.2.0)
 pkgconfig              2.0.3      2019-09-22 [1] CRAN (R 4.2.0)
 plotly                 4.10.0     2021-10-09 [1] CRAN (R 4.2.0)
 plyr                   1.8.7      2022-03-24 [1] CRAN (R 4.2.0)
 png                    0.1-7      2013-12-03 [1] CRAN (R 4.2.0)
 polyclip               1.10-0     2019-03-14 [1] CRAN (R 4.2.0)
 princurve            * 2.1.6      2021-01-18 [1] CRAN (R 4.2.0)
 prismatic              1.1.1      2022-08-15 [1] CRAN (R 4.2.0)
 progressr              0.10.1     2022-06-03 [1] CRAN (R 4.2.0)
 promises               1.2.0.1    2021-02-11 [1] CRAN (R 4.2.0)
 ps                     1.7.1      2022-06-18 [1] CRAN (R 4.2.0)
 purrr                  0.3.4      2020-04-17 [1] CRAN (R 4.2.0)
 R6                     2.5.1      2021-08-19 [1] CRAN (R 4.2.0)
 RANN                   2.6.1      2019-01-08 [1] CRAN (R 4.2.0)
 RColorBrewer           1.1-3      2022-04-03 [1] CRAN (R 4.2.0)
 Rcpp                   1.0.9      2022-07-08 [1] CRAN (R 4.2.0)
 RcppAnnoy              0.0.19     2021-07-30 [1] CRAN (R 4.2.0)
 RcppEigen              0.3.3.9.2  2022-04-08 [1] CRAN (R 4.2.0)
 RCurl                  1.98-1.8   2022-07-30 [1] CRAN (R 4.2.0)
 rematch2               2.1.2      2020-05-01 [1] CRAN (R 4.2.0)
 reshape2               1.4.4      2020-04-09 [1] CRAN (R 4.2.0)
 reticulate           * 1.28       2023-01-27 [1] CRAN (R 4.2.0)
 rgeos                  0.5-9      2021-12-15 [1] CRAN (R 4.2.0)
 rjson                  0.2.21     2022-01-09 [1] CRAN (R 4.2.0)
 rlang                  1.1.1      2023-04-28 [1] CRAN (R 4.2.0)
 rmarkdown              2.16       2022-08-24 [1] CRAN (R 4.2.0)
 ROCR                   1.0-11     2020-05-02 [1] CRAN (R 4.2.0)
 rpart                  4.1.16     2022-01-24 [1] CRAN (R 4.2.1)
 rprojroot              2.0.3      2022-04-02 [1] CRAN (R 4.2.0)
 rstudioapi             0.14       2022-08-22 [1] CRAN (R 4.2.0)
 Rtsne                  0.16       2022-04-17 [1] CRAN (R 4.2.0)
 rvest                  1.0.3      2022-08-19 [1] CRAN (R 4.2.0)
 S4Vectors            * 0.34.0     2022-04-26 [1] Bioconductor
 sandwich               3.0-2      2022-06-15 [1] CRAN (R 4.2.0)
 scales                 1.2.1      2022-08-20 [1] CRAN (R 4.2.0)
 scattermore            0.8        2022-02-14 [1] CRAN (R 4.2.0)
 scLANE               * 0.7.3      2023-09-20 [1] Bioconductor
 sctransform            0.3.4      2022-08-20 [1] CRAN (R 4.2.0)
 sessioninfo            1.2.2      2021-12-06 [1] CRAN (R 4.2.0)
 Seurat               * 4.1.1      2022-05-02 [1] CRAN (R 4.2.0)
 SeuratObject         * 4.1.0      2022-05-01 [1] CRAN (R 4.2.0)
 shape                  1.4.6      2021-05-19 [1] CRAN (R 4.2.0)
 shiny                  1.7.2      2022-07-19 [1] CRAN (R 4.2.0)
 SingleCellExperiment * 1.18.0     2022-04-26 [1] Bioconductor
 slingshot            * 2.4.0      2022-04-26 [1] Bioconductor
 sp                   * 1.5-0      2022-06-05 [1] CRAN (R 4.2.0)
 sparseMatrixStats      1.8.0      2022-04-26 [1] Bioconductor
 spatstat.core          2.4-4      2022-05-18 [1] CRAN (R 4.2.0)
 spatstat.data          3.0-0      2022-10-21 [1] CRAN (R 4.2.0)
 spatstat.geom          3.0-3      2022-10-25 [1] CRAN (R 4.2.0)
 spatstat.random        3.0-1      2022-11-03 [1] CRAN (R 4.2.0)
 spatstat.sparse        3.0-0      2022-10-21 [1] CRAN (R 4.2.0)
 spatstat.utils         3.0-1      2022-10-19 [1] CRAN (R 4.2.0)
 stringi                1.7.8      2022-07-11 [1] CRAN (R 4.2.0)
 stringr                1.4.1      2022-08-20 [1] CRAN (R 4.2.0)
 SummarizedExperiment * 1.26.1     2022-05-01 [1] Bioconductor
 survival               3.4-0      2022-08-09 [1] CRAN (R 4.2.0)
 svglite                2.1.0      2022-02-03 [1] CRAN (R 4.2.0)
 systemfonts            1.0.4      2022-02-11 [1] CRAN (R 4.2.0)
 tensor                 1.5        2012-05-05 [1] CRAN (R 4.2.0)
 TH.data                1.1-1      2022-04-26 [1] CRAN (R 4.2.0)
 tibble                 3.1.8      2022-07-22 [1] CRAN (R 4.2.0)
 tidyr                  1.2.0      2022-02-01 [1] CRAN (R 4.2.0)
 tidyselect             1.1.2      2022-02-21 [1] CRAN (R 4.2.0)
 TMB                    1.9.1      2022-08-16 [1] CRAN (R 4.2.0)
 TrajectoryUtils      * 1.4.0      2022-04-26 [1] Bioconductor
 utf8                   1.2.2      2021-07-24 [1] CRAN (R 4.2.0)
 uwot                   0.1.14     2022-08-22 [1] CRAN (R 4.2.0)
 vctrs                  0.6.3      2023-06-14 [1] CRAN (R 4.2.0)
 viridis                0.6.2      2021-10-13 [1] CRAN (R 4.2.0)
 viridisLite            0.4.1      2022-08-22 [1] CRAN (R 4.2.0)
 webshot                0.5.3      2022-04-14 [1] CRAN (R 4.2.0)
 withr                  2.5.0      2022-03-03 [1] CRAN (R 4.2.0)
 xfun                   0.32       2022-08-10 [1] CRAN (R 4.2.0)
 xml2                   1.3.3      2021-11-30 [1] CRAN (R 4.2.0)
 xtable                 1.8-4      2019-04-21 [1] CRAN (R 4.2.0)
 XVector                0.36.0     2022-04-26 [1] Bioconductor
 yaml                   2.3.5      2022-02-21 [1] CRAN (R 4.2.0)
 zlibbioc               1.42.0     2022-04-26 [1] Bioconductor
 zoo                    1.8-10     2022-04-15 [1] CRAN (R 4.2.0)

 [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /Users/jack/Desktop/Python/science/venv/bin/python
 libpython:      /usr/local/opt/python@3.8/Frameworks/Python.framework/Versions/3.8/lib/python3.8/config-3.8-darwin/libpython3.8.dylib
 pythonhome:     /Users/jack/Desktop/Python/science/venv:/Users/jack/Desktop/Python/science/venv
 virtualenv:     /Users/jack/Desktop/Python/science/venv/bin/activate_this.py
 version:        3.8.16 (default, Dec  7 2022, 01:36:11)  [Clang 14.0.0 (clang-1400.0.29.202)]
 numpy:          /Users/jack/Desktop/Python/science/venv/lib/python3.8/site-packages/numpy
 numpy_version:  1.23.5
 
 NOTE: Python version was forced by use_python function

──────────────────────────────────────────────────────────────────────────────